#include "iostream"
#include "cmath"
#include "ctime"
#include "random"

#include "../lib/lbfgs.h"

using std::clog;
using std::endl;

#define M 120
#define N 100

inline lbfgsfloatval_t scalernd(lbfgsfloatval_t l, lbfgsfloatval_t h)
{
    return (h-l)*rand()*1.0/RAND_MAX + l;
}

inline void vecrnd(lbfgsfloatval_t *a, lbfgsfloatval_t l, lbfgsfloatval_t h, int n)
{
    srand(time(nullptr));
	for (size_t i = 0; i < n; i++)
	{
		a[i] = (h-l)*rand()*1.0/RAND_MAX + l;
	}
    return;
}

inline lbfgsfloatval_t max_diff(const lbfgsfloatval_t *a, const lbfgsfloatval_t *b, int n)
{
	lbfgsfloatval_t tmp, max = -1.0;
	for (int i = 0; i < n; i++)
	{
        tmp = sqrt((a[i] - b[i])*(a[i] - b[i]));
		max = tmp>max?tmp:max;
	}
	return max;
}

class TEST_FUNC
{
public:
	TEST_FUNC();
	virtual ~TEST_FUNC();

	static lbfgsfloatval_t _Func(void *instance, const lbfgsfloatval_t *x, lbfgsfloatval_t *g,
		const int n, const lbfgsfloatval_t step)
	{
		return reinterpret_cast<TEST_FUNC*>(instance)->Func(x, g, n, step);
	}

	lbfgsfloatval_t Func(const lbfgsfloatval_t *x, lbfgsfloatval_t *g,
		const int n, const lbfgsfloatval_t step);

	static int _Progress(void *instance, const lbfgsfloatval_t *x, const lbfgsfloatval_t *g, const lbfgsfloatval_t fx,
		const lbfgsfloatval_t xnorm, const lbfgsfloatval_t gnorm, const lbfgsfloatval_t step, const lbfgs_parameter_t param,
		int n, int k, int ls)
	{
		return reinterpret_cast<TEST_FUNC*>(instance)->Progress(x, g, fx, xnorm, gnorm, step, param, n, k, ls);
	}

	int Progress(const lbfgsfloatval_t *x, const lbfgsfloatval_t *g, const lbfgsfloatval_t fx,
		const lbfgsfloatval_t xnorm, const lbfgsfloatval_t gnorm, const lbfgsfloatval_t step, const lbfgs_parameter_t param,
		int n, int k, int ls);

    static void _Precondition(void *instance, const lbfgsfloatval_t *x, const lbfgsfloatval_t *g, const lbfgsfloatval_t *d,
        lbfgsfloatval_t *d_pre, const int n)
	{
		return reinterpret_cast<TEST_FUNC*>(instance)->Precondition(x, g, d, d_pre, n);
	}

	void Precondition(const lbfgsfloatval_t *x, const lbfgsfloatval_t *g, const lbfgsfloatval_t *d,
        lbfgsfloatval_t *d_pre, const int n);

	int Routine();
private:
    lbfgsfloatval_t **kernel;
	lbfgsfloatval_t *x, *m, *mid, *d;
};

TEST_FUNC::TEST_FUNC()
{
	kernel = new lbfgsfloatval_t* [M];
    for (size_t i = 0; i < M; i++)
    {
        kernel[i] =  new lbfgsfloatval_t [N];
    }
    
    x = new lbfgsfloatval_t [N];
    m = new lbfgsfloatval_t [N];
    mid = new lbfgsfloatval_t [M];
    d = new lbfgsfloatval_t [M];

    srand(time(nullptr));
    for (size_t i = 0; i < M; i++)
    {
        for (size_t j = 0; j < N; j++)
        {
            kernel[i][j] = scalernd(-2.0, 2.0);
        }
    }

    vecrnd(m, 1.0, 2.0, N);

    for (size_t i = 0; i < M; i++)
    {
        d[i] = 0.0;
        for (size_t j = 0; j < N; j++)
        {
           d[i] += kernel[i][j] * m[j];
        }
    }
}

TEST_FUNC::~TEST_FUNC()
{
    for (size_t i = 0; i < N; i++)
    {
        delete[] kernel[i];
    }
    delete[] kernel;

    delete[] x;
    delete[] m;
    delete[] mid;
    delete[] d;
}

lbfgsfloatval_t TEST_FUNC::Func(const lbfgsfloatval_t *x, lbfgsfloatval_t *g,
		const int n, const lbfgsfloatval_t step)
{
    for (size_t i = 0; i < M; i++)
    {
        mid[i] = 0.0;
        for (size_t j = 0; j < N; j++)
        {
            mid[i] += kernel[i][j] * x[j];
        }
        mid[i] -= d[i];
    }

    lbfgsfloatval_t fx = 0;
    for (size_t i = 0; i < N; i++)
    {
        fx += mid[i]*mid[i];
    }
    
    for (size_t j = 0; j < N; j++)
    {
        g[j] = 0.0;
        for (size_t i = 0; i < M; i++)
        {
            g[j] += kernel[i][j]*mid[i];
        }
    }
	return fx;
}

int TEST_FUNC::Progress(const lbfgsfloatval_t *x, const lbfgsfloatval_t *g, const lbfgsfloatval_t fx,
		const lbfgsfloatval_t xnorm, const lbfgsfloatval_t gnorm, const lbfgsfloatval_t step, const lbfgs_parameter_t param,
		int n, int k, int ls)
{
	clog << "\riteration times: " << k << " fx = " << fx << " gnorm/xnorm = " << gnorm/(xnorm>1.0?xnorm:1.0);
	return 0;
}

void TEST_FUNC::Precondition(const lbfgsfloatval_t *x, const lbfgsfloatval_t *g, const lbfgsfloatval_t *d,
    lbfgsfloatval_t *d_pre, const int n)
{
    for (size_t j = 0; j < N; j++)
    {
        d_pre[j] = 0.0;
        for (size_t i = 0; i < M; i++)
        {
            d_pre[j] += kernel[i][j]*kernel[i][j];
        }
        d_pre[j] = d[j]/d_pre[j];
    }
    return;
}

int TEST_FUNC::Routine()
{
	lbfgsfloatval_t fx;

	lbfgs_parameter_t self_para;
	lbfgs_parameter_init(&self_para);
	self_para.epsilon = 1e-7;
    self_para.linesearch = LBFGS_LINESEARCH_BACKTRACKING_ARMIJO_QUAD;
    self_para.ftol = 1e-3;

	int ret = lbfgs(N, x, &fx, _Func, _Progress, this, &self_para, _Precondition);

	clog << endl << "L-BFGS optimization terminated with status: " << lbfgs_strerror(ret) << endl;
    clog << "Maximal difference = " << max_diff(x, m, N) << endl;
	return ret;
}

int main(int argc, char const *argv[])
{
	TEST_FUNC test;
	test.Routine();
	return 0;
}